## Replicate MVN simulations
## Last change: 6 April 2012
## Last used: 6 June 2012

library(MCMCpack)
library(RItools)
library(car)

## load core functions:
source("code/mahal.R")
source("code/seqblock1.R")
source("code/seqblock2k.R")
## create datasets
source("code/simmer.R")
set.seed(926)
simmer(times=100, outfileDir="data/simsMVN/r0/")
simmer(times=100, r=.8, outfileDir="data/simsMVN/r8/")
simmer(times=100, r=.6, whenoutlier = 2, outfileDir="data/simsMVN/r6outlier2/")
simmer(times=100, r=.6, whenoutlier = 20, outfileDir="data/simsMVN/r6outlier20/")
simmer(times=100, r=.6, whenoutlier = 35, outfileDir="data/simsMVN/r6outlier35/")
simmer(times=100, r=.8, covs = "contbim", outfileDir="data/simsMVN/r8bimodal/")

source("code/calcBalancePrecision.R")
source("code/calcBPdirectory.R")

## data.frames of all balance and precision statistics:
set.seed(817)
balprec.r0 <- calcBPdirectory("data/simsMVN/r0/", comp.rands = 100)
balprec.r8 <- calcBPdirectory("data/simsMVN/r8/", comp.rands = 100)
balprec.r6outlier2 <- calcBPdirectory("data/simsMVN/r6outlier2/", comp.rands = 100)
balprec.r6outlier20 <- calcBPdirectory("data/simsMVN/r6outlier20/", comp.rands = 100)
balprec.r6outlier35 <- calcBPdirectory("data/simsMVN/r6outlier35/", comp.rands = 100)
balprec.r8bimodal <- calcBPdirectory("data/simsMVN/r8bimodal/", comp.rands = 100)

save(balprec.r0, balprec.r8, balprec.r6outlier2, balprec.r6outlier20, balprec.r6outlier35, balprec.r8bimodal, file = "data/balprecAllMVN.RData")

load("data/balprecAllMVN.RData")
